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We discuss some mathematical aspects of the problem of inverting gravitational field data to 
extract the underlying mass distribution. While the forward problem of computing the gravity field 
from a given mass distribution is mathematically straightforward, the inverse of this forward map 
has some interesting features that make inversion a difficult problem. In particular, the forward 
map has an infinite-dimensional kernel which makes the inversion fundamentally non-unique. We 
characterize completely the kernels of two gravitational forward maps, one mapping mass density to 
the Newtonian scalar potential, and the other mapping mass density to the gravity gradient tensor, 
which is the quantity most commonly measured in field observations. In addition, we present some 
results on unique inversion under constrained conditions, and comment on the roles the kernel of the 
forward map and non-uniqueness play in discretized approaches to the continuum inverse problem. 

PACS numbers: 03.67.-a, 03.65.Ud, 04.70.Dy, 04.62. +v 



Weighing the shape of a gravitating body 

More than thirty years ago, Mark Kac asked "Can you hear the shape of a drum?" meaning: do two distinct 
planar domains always have distinct spectra of eigenvalues for their respective Laplace operators (acting on 
functions) with the usual Dirichlet (or Neumann) boundary conditions? If the answer is yes, the shape of a 
"drum" can be inferred by hearing its spectrum (characteristic sound), if the answer is no, then two distinctly 
shaped drums may have identical spectra (in which case they are called "isospectral domains" ) [2, HI • 

Kac's article [l[ stimulated a long line of research which eventually settled his question in the negative: There do 
exist isospectral domains (and, more generally, isospectral Riemann surfaces and isospectral Riemannian manifolds 
in higher dimensions) which are not isometric. In other words, the spectral inverse problem is ill-defined, subject 
to a fundamental ambiguity which can be precisely characterized Q. 

A similar ambiguity plagues the gravitational inverse problem, that is, the problem of inferring the precise 
shape of a mass distribution by observing its distant gravitational field. 

The gravitational inverse problem is the problem of inverting the gravitational forward map, which we take 
to be a map sending a compact supported mass distribution to a gravity observable: in practice, the observable 
could be either the Newtonian gravitational potential or gravity gradients. 

More precisely, and focusing on the gravity potential $ for the moment, what we will mean by the gravitational 
inverse problem is the following: Given a spherical region Br = {f : \r\ < R} of radius R in K 3 , and a solution 
$(r) (the gravitational potential in free space) of the Laplace equation V 2 <I> = outside the region Br (i.e. for 
\f\ > R) which vanishes at infinity, find a mass density distribution p(f') supported inside Br which gives rise to 
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$(r) in the exterior region outside Br. In plainer language, find a p(f') with support inside Br such that 

$(r) = -G [ f (r ? rfV for r>i?. (1) 



Kernel of the forward map onto the gravitational potential 

Equation (1) of course represents the unique solution to the "forward problem" where one searches for a solution 
<f> to V 2( f> = AirGp with vanishing boundary conditions at infinity. Of key interest is the "kernel" of this (linear) 
forward map, i.e. the set of mass distributions p supported inside Br that are mapped to a potential <f> via Eq. (1) 
which identically vanishes outside Br. 

Theorem 1 : The kernel of the forward map Eq. (1) mapping mass distributions p supported in Br to solutions 
of Laplace equation outside the region Br (i.e. for \r\ > R) is precisely functions p satisfying 

P = V 2 X , (2) 

where x(r) is any (sufficiently smooth) function on R 3 with support inside Br (i.e. \(^) = for r > i?). In other 
words, if p is a solution of the inverse problem for a given exterior potential <f>, then p + V 2 x is also a solution for 
any \ £ Co" (Br), where a is a sufficiently large integer. Normally, a > 2 should be sufficient, but smoothness is 
not a key issue; in particular, \ can even be a distribution if point-mass (delta-function) singularities need to be 
allowed in the problem. 

Proof in one direction is easy: Every function in the kernel is given by the forward image of a function of the 
kind Eq. (2). To prove this, let $ be a function belonging to the kernel, i.e. let $ vanish outside Br. Put 




Then \ € Co* (Br) and p = V 2 \ satisfies the Laplace equation V 2 <& = 4nGp everywhere (with vanishing boundary 
conditions at infinity). Therefore, $ satisfies Eq. (I) with this p, which is what we needed to prove. 

Conversely, let p be a density distribution supported inside Br such that p — V 2 x for some x G Co" (Br). 
Then, according to Eq. (I), the gravitational potential $ which the forward map sends p onto satisfies 

$(r) = - G [ ^y^J d 3 r' for r > R . (3) 
Jb r \r-r'\ 

To show that the right hand side of Eq. (3) is in the kernel of the forward map, i.e., that it vanishes for r > R, 
use Green's identity: 

l B V*v-v*u)*T-f a (u%;-vg)*, <«) 

where B is any region bounded by the surface dB, and U, V are arbitrary functions on R 3 . Applying Eq. (4) 
with B taken as the region Br, U(r') = l/\r — r'\, and V(r') = x(0; an d noting that V 2 (I/|r — r'\) — when 
r > R and r' < R, it immediately follows that the right hand side of Eq. (3) vanishes outside Br (i.e. for r > R). 
This completes the proof of Theorem I . 



A geometric interpretation of the kernel: 

One way to conceptualize the kernel of the gravitational (potential) forward map is to note that the geometric 
freedom of choice in the inverse "datum" $(r) is that of choosing an arbitrary function on the two-sphere S 2 : 

Theorem 2 : Let $(r) be any solution of the free-space Laplace equation V 2 $ = outside the region Br which 
vanishes at infinity, as in the formulation of the gravitational inverse problem. Then $ is completely determined 
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outside Br by its values on any two-sphere Sri. of radius R\ > R (or, more generally, by its values on any closed 
surface which encloses Br). 

Proof: This is really a restatement of a standard result in potential theory (uniqueness of solutions to the 
Dirichlet problem): There exists a unique Green's function G(r, fo), defined for r, r outside Br 11 such that 
G satisfies V 2 G(r, f ) — S(r — fo) and vanishes for f G Sr 1 and for f G Soo [as discussed, e.g., in [4(, for the 
two-sphere Sr ± G(f,fo) can be constructed explicitly using the classic "method of images"]. Plugging such a 
G(r,f ) into Eq. (4) as V and taking <I>(r) as U and the region B as the region outside Br 1 we obtain, by virtue 
of the vanishing boundary conditions at infinity, 

*(*>)=/ m^^da. (5) 

Therefore $ everywhere outside Br 1 is determined uniquely by its values on the two-sphere . 

We can now understand the kernel Eq. (2) in the following way: Since the data for $ consist of the values of 
a function defined on a two-surface Sr 1 , we can infer from these data uniquely at the most another function of 
two variables, and not the full three-dimensional density field p{f). In fact, the forward kernel (or the ambiguity 
in the corresponding inversion) as described by Eq. (2) corresponds precisely to this geometric statement. 



Weighing the shape of a body of known radial density 

One might hope that practical (physical) prior constraints on the three-dimensional density distribution p(f) 
might make it uniquely recoverable from its far-zone gravity field despite the fundamental non-uniqueness of 
the inverse problem. For example, we want the density to be positive everywhere, which is a requirement that 
constrains the ambiguity Eq. (2) to some extent. However, simple spherically-symmetric counterexamples show 
that positivity is not a sufficiently strong constraint to help provide us with a unique inversion. As the next step 
in a series of physically-reasonable constraints on p, we might assume a known positive radial density distribution 
with profile p(f) = po(r) > distributed on some arbitrary compact three-dimensional region D in R 3 . Put 
another way, such a density profile represents a body of arbitrary shape carved out of a spherically symmetric 
(hence spherical) mass distribution. Again, counterexamples based on hollow spherical shells show that this is not 
quite enough for unique inversion. Nevertheless, it turns out that if we further constrain the region D such that 
it is connected and has no "holes" (i.e., if D is topologically a ball), and, furthermore, if D is "radially convex" 
in a sense made precise below, then unique inversion is possible: 

Theorem 3 : Let D be compact region in M 3 such that its boundary dD is a connected and simply-connected 
surface (in other words, dD is a topological two-sphere) which is radially convex in the sense that any straight 
line in R 3 passing through the center-of-mass of the volume D intersects dD at precisely two points. Assume 
that D is filled with material of a known non-negative mass density p{f) which, when it is nonzero, is distributed 
spherically-symmetrically with respect to the coordinate origin given by the center of mass. That is, if r lies inside 
D, then p{f) = po(r) > 0, and if r is outside D, then p(f) — 0. Under these conditions, D itself (or, equivalently, 
its boundary dD) is uniquely recoverable from the far-zone gravity field of this radial density distribution. ^ 

This result is not too surprising in view of Theorem 2, since the specification of dD entails just a single real 
function on the two-sphere S 2 (measuring just how much we need to deform S 2 in order to stretch it onto dD). 
The forward map Eq. (1) can then be interpreted as a nonlinear map from real functions on S* 2 (representing the 
deformations of S 2 needed to obtain dD) to real functions on S 2 (representing the values of the potential <3? on 
Sr), and we will now show that this map is locally one-to-one. 

Proof of Theorem 3: The main idea of the proof is simple: explicitly write down, in terms of spherical-harmonic 



t The assumption that dD is a topological two-sphere is redundant since it follows from the assumption of radial convexity as 
formulated in the theorem. However, it is perhaps useful to emphasize this assumption in a redundant statement since the theorem 
is certainly false without it. 



4 



coefficients, the forward transform mapping the "shape function" of dD to the exterior potential $, and show 
that the derivative of this nonlinear forward map is nonsingular. The result then follows from the inverse function 
theorem as generalized to infinite-dimensional spaces [5]. In this paper we will give a detailed proof that the 
forward map has nonsingular derivative at the point (shape) which corresponds to a perfect sphere, so the result 
holds for shapes which are nearby distortions of a perfect sphere (in other words, we will explicitly prove that 
the forward map is invertible in some open neighborhood of the perfect sphere in the space of all shapes D which 
satisfy the conditions of the theorem). This case covers most planetary bodies at the levels of resolution we 
are interested in. Nevertheless, the statement that the forward map is nonsingular everywhere remains valid, 
although we are not going to give an explicit proof of it here. The proof of this more general case is substantially 
similar apart from requiring more careful estimates. 

To proceed with the proof, introduce a spherical coordinate system (r, 9, <p) centered at the center of mass of 
the volume D. In this coordinate system, let n(8, </>) denote the half-line which starts at the origin and expands 
outward in the direction (6,4>). Let ip(6,(j)) be the (positive) function which gives the length of the radial vector 
which starts at the origin and ends at the intersection point of the line n(8, <f) with the boundary dD as (9, cj>) 
ranges over the unit two-sphere S* 2 of all possible directions [by the radial convexity assumption, there exists a 
unique such intersection point for each direction (9, </))]. The function ^){9, </>) can then be taken to be the "shape 
function" which specifies D, and, explicitly, we can write 

D = {(r,(l) \r< V(O)}, (6) 

and 



dD = {(r,Sl) |r = ^(n)}, (7) 

where we introduced the short- hand notation il = (6,<p) for the angular coordinates. The exterior gravitational 
potential $ can be expanded in spherical harmonics [4]: 

$(r> - dl ™^0T- iovr>R> max , (8) 

where we can regard {d/ m } = D as an infinite sequence (vector) of "observables" which completely describes 
the data for the inverse problem in view of Theorem 2. On the other hand, according to Eq. (1), for r > R > 
maxn ip(Q) we have 

r d 3 r' 
m = G / Poir')^^ 
J r ><-4>{n<) \r-r'\ 

r*l>(. a ' 



= G po(r')r ,2 dr' £ Y^Y^n') dV 

J o Js* 



GE/-N%P> (9) 



I, in 



where 



and 



Hn(w)= / Po (r)r n dr, (10a) 



/ Im M = / Y; m (n) wo)) dn (10) 
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is a vector functional of the shape function i/j which represents the forward map in the same way as D = {d; m } 
represents the data. In fact, introducing the notation F[i/j] = {/; m [V']} an d combining Eqs. (8) and (9), the 
forward equation for the shape function if) takes the simple form 

m = ^ d . (ii) 

(Due to our choice of origin as the center of mass, both fi m and di m vanish for 1 = 1, but this fact will not be 
of any consequence in what follows.) It is also convenient to introduce a coordinatization of the space of shape 
functions via a spherical harmonic expansion 

lP(Q.) = J2simYlm(fy , (12) 

and consider the coordinate vector S = {s; m } as the representation of the function ip(Q,). In this coordinate 
system the forward map Eq. (11) takes the form 

F[S] = i D , (13) 

where 

/ Im [s] = r^(o) Mi+2 (^2 s pM^ dn ■ ( 14 ) 

Assume now, contrary to the conclusion of Theorem 3, that two distinct domains D\ and D% constrained as in the 
statement of the theorem give rise to identical external gravitational potentials when filled with the given radial 
density distribution po(r). First of all, since the monopole and dipole moments of the two mass distributions 
must agree, they must have the same center of mass, therefore we can set up a common spherical coordinate 
system for both volumes with their shared center of mass chosen as the origin of coordinates. It then follows that 
there exist two distinct shape functions if)\ and ip2, corresponding to the two distinct volumes D\ and D2, which 
satisfy Eq. (11) with the same data D; in other words 

F[H = F[H . (15) 

We will now show that Eq. (15) is impossible as long as ipi and ip2 belong to some fixed open neighborhood of 
a perfect sphere {ip(Q) = &o = const} in the infinite-dimensional nonlinear function space of all ^'s. Using the 
inverse function theorem as generalized to such infinite-dimensional manifolds Q , it is sufficient to show that the 
derivative of the map F at the point tj)(Cl) = ao is a nonsingular linear map. In general, at an arbitrary point 
ip = ipa(£l), this derivative is given by 

(F'[V>o(0)] • tyOim = / Ybnffl Po (Vo(O)) M^Y +2 S^(Cl)dn , (16) 
Js 2 

where F'[^o(^)] denotes the derivative evaluated at the point ip = ipo, acting (as a linear map) on the tangent 
vector (linear perturbation) Sip, and we have used Eq. (10) to derive this explicit form. Specializing to the perfect 
sphere ipo(Q) = clq(= const) and using the coordinate representation [cf. Eq. (12)] 

ty(fi) = KMty , SS = {Ss pq } , (17) 

Eq. (16) takes the form 



(F'[a ] • SS) lm = p (a )a l+2 / Y^Q) 



s- 



Lp, q 

1+2 



po(ao) a Ssim , (18) 
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where we made use of the fact that the Y/ m 's form an orthonormal basis for L 2 (S 2 ). According to Eq. (18), the 
derivative F'[a ] is a diagonal linear map with only nonzero entries (eigenvalues) on the diagonal; therefore, F'[a ] 
is clearly nonsingular. This completes the proof of Theorem 3. 



The kernel of the forward map onto gravity gradient observables 

The gravitational gradient tensor is (apart from a minus sign) simply the (symmetric) tensor of second deriva- 
tives of the potential $ in a cartesian coordinate system: 

<9 2 $ 

TlJ = 'dx^dxl ■ (19) 

So, for example, we have 

<9 2 $ d 2 $ 

Txx = ~dx^ ' Tyz = ~^y~d~z (20) 

etc. Independently of coordinates, the gradient tensor can be defined as the double covariant derivative VV$ (in 
general relativity, corresponds to the Riemann curvature tensor Rowj describing tidal gravitational forces). 
One can also define the gradient tensor explicitly in terms of the source mass distribution as: 

T«(f) = G ( PMW-vW-g-^-vr] d s y > (21) 

J \ x y\ 

where all coordinates are cartesian. The gradient tensor is a particularly useful observable in precision gravimetry 
since it is better isolated from local non-gravitational acceleration noise compared to other observables, and a 
large roster of instruments (gradiometers) are available for measuring it. 

In practical applications, one often works with a coordinate system where z is the vertical coordinate pointing 
up from the Earth's center, and the observable of interest is the x — y projection of the gradient tensor in an 
infinitesimally small neighborhood (tangent plane to the Earth's spherical surface) around x = y = 0: 



T T 

-l xx -l xy 

T T 

J-xy J-yy 



(22) 



Typically, a gradiometer takes two kinds of measurements: the component M x = 2T xy ("crossline"), and the 
combination M + = T xx — T yy (inline). The choice of cartesian x, y coordinates is arbitrary upto a rotation R, 
and T transforms under rotations as 

T — > RTR 1 . (23) 

Neither T nor its crossline or inline components are invariant under rotations, but Tr(T) and Det(T) are invariants. 
In particular, the Euclidean norm of (M + ,M X ) 



' Ml + M 2 . = v/Tr(T) 2 -4Dct(T) (24) 
is an invariant. More specifically, it is easy to compute that (M + ,M X ) transforms under a rotation 

/cos^ -sin<A 

\sva.v cos 9 J v ' 



according to the rule 



mJ ^ R20 [m 



; ^U + l- (26) 
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Now, define a new observable 

V = Ml + Ml . (27) 

In other words, V is defined at every point of M 3 by setting up local coordinates x, y, z such that the z-axis goes 
through the origin, computing the gradient tensor T as in Eqs. (19) and (22) at that point, and then calculating 
the norm square of the observable (M + ,M X ). The key features of V are that (i) it is invariant under rotations, 
therefore uniquely defined independently of the choice of coordinates, and (ii) it obeys the following lemma: 

Lemma: V vanishes at a point if and only if (-M+ , M x ) vanishes there for any allowed choice of coordinates 
x,y, z. 

Proof: By Eq. (11), V is invariant and equal to the norm-square of (M + ,M X ) for any choice of coordinates. 

Therefore the kernel of the coordinate-dependent observable (M+, M x ) is precisely the kernel of the nonlinear 
but coordinate-independent observable V . 

Theorem 4: Let be any analytic functional (linear or nonlinear) of the gravity potential $ (such as V), 

and let S be any analytic 2-surface lying outside the spherical region Br = {f: \r\ < R} in R 3 (such as a sphere 
of radius > R). Then, if F[<&\ vanishes in any open (two-dimensional) neighborhood on S, it vanishes identicallly 
on all of S. 

Proof: This just follows from analyticity: $ is an analytic function outside the spherical region Br = {r : \r\ < R} 
in K 3 since it satisfies the homogeneous Laplace equation there (analyticity follows from standard elliptic regularity 
theorems @). Therefore, F[$] is analytic there and so is its restriction to S since S is analytic. Thus vanishing 
on any open subset is equivalent to vanishing identically on S. 

Corollary: If V vanishes in any two-dimensional patch, no matter how small, on any analytic observation surface 
S lying outside the spherical region Br = {r : \r\ < R}, then it vanishes identicallly on all of S. 

This corollary further illustrates the fact that the kernel of the gravity-gradient observable (M+ , M x ) is precisely 
the kernel of the observable V, not only locally but also globally. 

Theorem 5: Let S be a sphere of radius > R. Then V vanishes on S if and only if $ is spherically symmetric 
(a function of the radius r only), and hence V vanishes identically everywhere outside the spherical region 
B R = {r : \r\ < R} in M 3 . 

Proof: The if part is a simple calculation: It is straightforward to compute that both M+ and M x vanish for 
a radial (monopole) potential function $(?•). For the converse, it is easy to see that if V vanishes on a sphere S 
then this implies that $ is constant on S. But this implies, according to Eq. (5) (or just by the uniqueness of 
solutions to the Dirichlet problem), that $ is a radial function of monopole type: 

$(r) = - , (28) 
r 

where C is a constant. 

We can now completely characterize the kernel of the forward map from the mass density to the gravity gradient 
observables (M + , M x ) : 

Theorem 6 : The kernel of the forward map mapping mass distributions p supported in Br to gravity gradient 
observables (M + ,M X ) outside the region Br (i.e. for \r\ > R) is precisely functions p supported in Br satisfying 

p = Po + V 2 X , (29) 

where x(r) is any (sufficiently smooth) function on R 3 , and po is any spherically symmetric function, both 
supported inside Br (i.e., both x(f) and po(r) vanish for r > R). 

Proof: By Theorem 5, the kernel of the map from mass distributions to outside gradients (M + ,M X ) consists 
of those mass distributions that give rise to spherically symmetric potentials <& outside Br. If $ is spherically 
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symmetric outside Br, then consider $s, the spherical average of $ (i.e. <&s(%) = the average of $ on the sphere 
centered at and passing through x). Then $5 is spherically symmetric everywhere and coincides with $ outside 
Br. Therefore x = $ — $5 vanishes outside Br, and thus 

V 2 $ = V 2 x + V 2 $ s . (30) 

Since $s is everywhere spherically symmetric, so is V 2 $s, and Theorem 6 follows. 



Kernel of the gravitational forward map and discretization 

In practice, gravity inversion is a discrete problem because (i) the measurements of <E> (or of the gradients) are 
finitely many and discretely distributed in space, and (ii) more problematically, the model for the mass distribution 
p is some discretized approximation to a continuous distribution. In all contexts, the characterization of p would 
be a finite list of parameters which uniquely specify p in some generally non-linear fashion. For example, these 
p parameters could be the masses, locations, and shape parameters of a finite number of tectonic plates in a 
geophysical model of the Earth's crust. Or, as discussed above, they could be masses of finite blocks into which 
we divide the source distribution in discretizing it. Or, more straightforwardly, they could be the masses of TV 
point-mass centers distributed throughout the source region, approximating with a discrete configuration the true 
mass distribution in the limit N — > 00. 

In general, the practical, the discretized gravitational inverse problem is the problem of inverting some (generally 
nonlinear) forward map: 

F: fo} _►{$,}, Z^Ffa] , (31) 

where pj are finitely many parameters specifying the mass distribution, and $i are the measurements. It would be 
conceptually salutary to have the fundamental non-uniqueness in the gravitational inverse problem (the kernel of 
the froward map) described by Theorem 1 to fall out of the formulation Eq. (31 ) in a natural way. For example, 
when we simulate a slab of soil using some large number of mass centers regularly placed at fixed lattice points 
inside the slab, the parameters pj are simply the point masses rrij assigned to each center at lattice location j. 
In this case, the forward map F is in fact linear: 

where the matrix F^ is the Green's function in Eq. (1) in discretized form: 

F " - - G inhri ■ (33) 

with f\ being the locations where the measurements $i = &(fl) are collected. Consider first, for simplicity, a 
scenario in which we are sampling $ at the same number of points N as the number of mass centers in the 
discretization. In other words, F is now a square N x N matrix. In view of Theorem 1 characterizing the kernel 
of the forward map, one might expect F to be singular, with the null space corresponding to a discretized version 
of the kernel, i.e., a discrete approximation to functions of the form V 2 x with x supported inside the slab. It 
turns out, however, that the matrix F given by Eq. (33) is in fact generically nonsingular. Moreover, the M x N 
matrix F with M measurement locations and N mass centers is also nonsingular, in the sense that generically it 
has maximal rank (i.e. trivial null space). 

It is in fact easy to see why this is so, because of the following result: 

Theorem 7 : Given a solution $(r) of the Laplace equation V 2 <1> = vanishing at infinity and defined for r > R, 
there exists at most one configuration {mj,fj} of finitely many point masses placed inside B R — {r < R} (i.e. 
with rrij e R and \fj\ < R) that can give rise to this $ for r > R. 

Proof: Suppose, on the contrary, that there are two configurations of point masses, {mj,fj} and {m' k ,fk'}, 
that produce the same $ for r > R. Subtract the second configuration from the first, and correspondingly 
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subtract the $ fields that they produce. Since the gravity field depends linearly on the mass distribution, what 
we obtain is a new configuration {mi, • • • ,mjv, -m'j, • • • — m' N ,,f{, • • • , r]v, f{' , • • • , rjv''} of point masses inside 
Br (unless there are some coincident point masses in the two collections, in which case one would simply subtract 
the corresponding masses and list the location only once), which produces a field $ that vanishes identically for 
r > R. Could this actually happen? It turns out the answer is no, unless $ is identically zero everywhere (and 
therefore the two original point- mass configurations are in fact identical). To see this, observe that <1> produced 
by a finite set of point masses is a real-analytic function in M 3 except at the locations of the point masses where 
it has singularities. Since $ vanishes for r > R and is analytic, it must vanish everywhere in R 3 except possibly 
at the mass centers. But if any of the mass centers had non-zero mass, we could choose points so close to that 
center that the contribution to $ from that center would overwhelm the contributions from any other centers 
(which are discretely spaced since there are finitely many). This clearly contradicts the fact that $ is identically 
zero in any small neighborhood of the chosen mass-center. Therefore, none of the mass centers can have nonzero 
mass; the two original configurations of point masses must be identical, and Theorem 7 is proved. 

Here is one way to understand the apparent conflict between Theorem 1 and Theorem 7: Consider the two 
spaces between which the forward map F acts: the space of density distributions p and the space of potentials 
$. Any discretization is an attempt to approximate these spaces via a sequence of finite-dimensional subspaces. 
For example, when we use N point masses, we have an iV-dimensional subspace of the space of all p, and as N 
gets larger and larger this subspace approximates the full space arbitrarily closely, in the sense that for any p , 
we can find a configuration of N point masses (with large enough N) which comes as close as we want to p (in 
some locally averaged sense). The same goes for the corresponding potentials $: given any solution $0, we can 
find potentials produced by N point masses that get arbitrarily close to $o as N — > oo. But the problem is that 
these approximating subspaces completely miss the kernel of the true forward map, which is the subspace of mass 
distributions (and corresponding potentials) given by {V 2 x | x G Co(-Br)}. The intersection of the approximating 
subspaces with this kernel subspace is the zero vector, for any finite N. This is (mathematically) the explanation 
for the apparent contradiction between Theorem 1 and Theorem 7. 

To resolve this apparent conceptual paradox, one might argue that we must choose the approximating finite 
dimensional subspaces in such a way that they fully intersect the kernel. But in practice, there is no feasible way to 
discretize the problem that makes sure this property holds. There is, however, a much simpler practical strategy 
out of this apparent paradox, and this is the strategy we advocate: Realize that true measurements in the real 
world always have instrumental noise. What this means is that two potentials are indistinguishable in practice if 
they differ everywhere by less than, say la worth (in some arbitrary units) of instrumental noise. Therefore, e.g. 
when we look for the intersection between the kernel and our discretized p-subspace with N point masses, what 
we are really looking for are all iV-point-mass configurations that produce a potential $ that differs from zero by 
less than la throughout the exterior region r > R. And in general there are many such configurations. We can see 
this, for example, in the matrix F^j of Eq. (33): in general this matrix turns out to be highly ill-conditioned (with 
very small determinant) with lots of eigenvalues close to zero, even though it has no exactly-zero eigenvalues. 
And the "approximately null" subspace spanned by the small-eigenvalued eigenspaces is precisely the discrete 
analogue of the kernel of the forward map; it is what corresponds to the subspace {V 2 x | X S Co(Br)} in this 
discretization. We expect a similar description for the discrete analogue of the forward map's kernel in any other 
practical discretization scenario. 
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